Load the necessary libraries
library(car) #for regression diagnostics
library(broom) #for tidy output
library(broom.mixed) #for tidy output
library(ggfortify) #for model diagnostics
library(sjPlot) #for outputs
library(knitr) #for kable
library(effects) #for partial effects plots
library(ggeffects) #for effects plots in ggplotjk
library(emmeans) #for estimating marginal means
library(MASS) #for glm.nb
library(MuMIn) #for AICc
library(tidyverse) #for data wrangling
library(DHARMa) #for assessing dispersion etc
library(glmmTMB) #for glmmTMB
library(performance) #for diagnostic plots
library(see) #for diagnostic plots
library(lme4) #for glmer
library(glmmTMB) #for glmmTMB
Some ornithologists were interested in the degree of sibling negotiations in owl chicks. Specifically, they wanted to explore how sibling negotiations were affected by feeding satiety and the sex of the parent returning to the nest. The ornithologists had accessed to a number of owl nests and were able to count (via recording equipment) the number of sibling negotiations (calls) that the owl chicks made when the parent returned to the nest.
We could hypothesise that the chicks might call more if they were hungry. As part of the investigation, the researchers were able to provided supplimentary food. As such, they were able to manipulate the conditions such that sometimes the chicks in a nest would be considered deprived of supplimentary food and at other times they were satiated.
As a parent returned, the researchers recorded the number of sibling negotiations (calls) along with the sex of the parent. Since the number of calls is likely to be a function of the number of chicks (the more chicks the more calls), the researchers also counted the number of siblings in the brood.
Each nest was measured on multiple occasions. Hence, we must include the nest as a random effect to account for the lack of independence between observations on the same set of siblings.
owls = read_csv('../public/data/owls.csv', trim_ws=TRUE)
##
## ── Column specification ───────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
## cols(
## Nest = col_character(),
## FoodTreatment = col_character(),
## SexParent = col_character(),
## ArrivalTime = col_double(),
## SiblingNegotiation = col_double(),
## BroodSize = col_double(),
## NegPerChick = col_double()
## )
glimpse(owls)
## Rows: 599
## Columns: 7
## $ Nest <chr> "AutavauxTV", "AutavauxTV", "AutavauxTV", "Autavau…
## $ FoodTreatment <chr> "Deprived", "Satiated", "Deprived", "Deprived", "D…
## $ SexParent <chr> "Male", "Male", "Male", "Male", "Male", "Male", "M…
## $ ArrivalTime <dbl> 22.25, 22.38, 22.53, 22.56, 22.61, 22.65, 22.76, 2…
## $ SiblingNegotiation <dbl> 4, 0, 2, 2, 2, 2, 18, 4, 18, 0, 0, 3, 0, 3, 3, 6, …
## $ BroodSize <dbl> 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5, 5,…
## $ NegPerChick <dbl> 0.8, 0.0, 0.4, 0.4, 0.4, 0.4, 3.6, 0.8, 3.6, 0.0, …
Let start by declaring the catogical variables and random effect as factors.
## Amount of Sibling negotiation (vocalizations when parents are absent)
## Foot treatment (deprived or satiated
## Sex of parent
## Arrival time of parent
## Nest as random
## Brood size offset
owls = owls %>% mutate(Nest =factor(Nest),
FoodTreatment = factor(FoodTreatment),
SexParent = factor(SexParent),
NCalls = SiblingNegotiation)
As the response represents counts (the number of calls), it would make sense to start by considering a Poisson model. We could attempt to model the response as the number of calls divided by the broodsize, but this would result in a response that has no natural distribution.
Instead, if we include broodsize as an offset, it will standardize the effects according to broodsize (similar to having divided the response by broodsize), yet retain the Poisson nature of the response.
The effects of offsets, unlike regular covariates, are not estimated. Rather they are assumed to be 1 (on the link scale). This means that since Poisson uses a log link, then the offset should be of a logged version of the broodsize.
Model formula: \[ y_i \sim{} \mathcal{Pois}(\lambda_i)\\ ln(\lambda_i) =\boldsymbol{\beta} \bf{X_i} + \boldsymbol{\gamma} \bf{Z_i} \]
where \(\boldsymbol{\beta}\) and \(\boldsymbol{\gamma}\) are vectors of the fixed and random effects parameters respectively and \(\bf{X}\) is the model matrix representing the overall intercept and effects of food treatment, sex of parent, arrival time (and various interactions) on the number of sibling negotiations. Brood size was also incorporated as an offset. \(\bf{Z}\) represents a cell means model matrix for the random intercepts associated with individual nests.
Perhaps we could start off by exploring the main fixed effects. To mimic the log-link, we will use a log-transformed y axis. Since there may well be zeros (no calls detected), we will use a psuedo log scale). We will also include the raw data (jittered and dodged)
ggplot(data = owls, aes(y = NCalls, x = FoodTreatment, color=SexParent)) +
geom_violin() +
geom_point(position=position_jitterdodge(jitter.height=0, dodge.width=1))+
scale_y_continuous(trans=scales::pseudo_log_trans())
Now, a similar plot separated for each nest.
ggplot(data=owls) +
geom_point(aes(y=NCalls, x=FoodTreatment, color=SexParent), position=position_dodge(0.5)) +
facet_wrap(~Nest)
It might also be worth establishing that there is a linear relationship between the number of calls and broodsize. Again, we will mimic the use of the log-link by transforming the axes.
ggplot(data = owls,aes(y = NCalls, x = BroodSize, color=SexParent)) +
geom_point() +
geom_smooth(method='lm') +
facet_grid(~FoodTreatment) +
scale_y_continuous(trans=scales::pseudo_log_trans()) +
scale_x_log10()
## `geom_smooth()` using formula 'y ~ x'
owls.glmer1 <- glmer(NCalls ~ FoodTreatment*SexParent + offset(log(BroodSize)) + (1|Nest), data=owls,
family=poisson(link='log'))
owls.glmer2 <- update(owls.glmer1, ~ . - FoodTreatment:SexParent)
AICc(owls.glmer1, owls.glmer2)
anova(owls.glmer1, owls.glmer2)
Conclusions:
owls.glmer1a <- owls.glmer1
owls.glmer1b <- update(owls.glmer1a, ~ . - (1|Nest) + (FoodTreatment|Nest))
owls.glmer1c <- update(owls.glmer1a, ~ . - (1|Nest) + (SexParent|Nest))
owls.glmer1d <- update(owls.glmer1a, ~ . - (1|Nest) + (FoodTreatment*SexParent|Nest))
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00681728 (tol = 0.002, component 1)
owls.allFit <- allFit(owls.glmer1d)
## Loading required namespace: dfoptim
## Loading required namespace: optimx
## bobyqa : [OK]
## Nelder_Mead :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 1.16117 (tol = 0.002, component 1)
## [OK]
## nlminbwrap : [OK]
## nmkbw :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00267498 (tol = 0.002, component 1)
## [OK]
## optimx.L-BFGS-B :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00238787 (tol = 0.002, component 1)
## [OK]
## nloptwrap.NLOPT_LN_NELDERMEAD : [OK]
## nloptwrap.NLOPT_LN_BOBYQA :
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.00642472 (tol = 0.002, component 1)
## [OK]
owls.allFit
## original model:
## NCalls ~ FoodTreatment + SexParent + (FoodTreatment * SexParent | Nest) + Foo...
## data: owls
## optimizers (7): bobyqa, Nelder_Mead, nlminbwrap, nmkbw, optimx.L-BFGS-B,nloptwrap.NLOPT_LN_N...
## differences in negative log-likelihoods:
## max= 0.739 ; std dev= 0.279
## Check which of the models are considered valid (OK)
is.OK <- sapply(owls.allFit, is, "merMod")
is.OK
## bobyqa Nelder_Mead
## TRUE TRUE
## nlminbwrap nmkbw
## TRUE TRUE
## optimx.L-BFGS-B nloptwrap.NLOPT_LN_NELDERMEAD
## TRUE TRUE
## nloptwrap.NLOPT_LN_BOBYQA
## TRUE
diff_optims.OK <- owls.allFit[is.OK]
lapply(diff_optims.OK,function(x) x@optinfo$conv$lme4$messages)
## $bobyqa
## NULL
##
## $Nelder_Mead
## [1] "Model failed to converge with max|grad| = 1.16117 (tol = 0.002, component 1)"
##
## $nlminbwrap
## NULL
##
## $nmkbw
## [1] "Model failed to converge with max|grad| = 0.00267498 (tol = 0.002, component 1)"
##
## $`optimx.L-BFGS-B`
## [1] "Model failed to converge with max|grad| = 0.00238787 (tol = 0.002, component 1)"
##
## $nloptwrap.NLOPT_LN_NELDERMEAD
## NULL
##
## $nloptwrap.NLOPT_LN_BOBYQA
## [1] "Model failed to converge with max|grad| = 0.00642472 (tol = 0.002, component 1)"
owls.glmer1d <- update(owls.glmer1d, control=glmerControl(optimizer='bobyqa'))
owls.glmer1c <- update(owls.glmer1c, control=glmerControl(optimizer='bobyqa'))
owls.glmer1a <- update(owls.glmer1a, control=glmerControl(optimizer='bobyqa'))
owls.glmer1b <- update(owls.glmer1b, control=glmerControl(optimizer='bobyqa'))
AICc(owls.glmer1a, owls.glmer1b, owls.glmer1c, owls.glmer1d)
owls.glmmTMB1 <- glmmTMB(NCalls ~ FoodTreatment*SexParent + offset(log(BroodSize)) + (1|Nest), data=owls,
family=poisson(link='log'), REML=FALSE)
owls.glmmTMB2 <- update(owls.glmmTMB1, ~ . - FoodTreatment:SexParent)
AICc(owls.glmmTMB1, owls.glmmTMB2)
anova(owls.glmmTMB1, owls.glmmTMB2)
Conclusions:
owls.glmmTMB1a <- update(owls.glmmTMB1, REML=TRUE)
owls.glmmTMB1b <- update(owls.glmmTMB1a, ~ . - (1|Nest) + (FoodTreatment|Nest))
owls.glmmTMB1c <- update(owls.glmmTMB1a, ~ . - (1|Nest) + (SexParent|Nest))
owls.glmmTMB1d <- update(owls.glmmTMB1a, ~ . - (1|Nest) + (FoodTreatment*SexParent|Nest))
AICc(owls.glmmTMB1a, owls.glmmTMB1b, owls.glmmTMB1c, owls.glmmTMB1d)
plot_model(owls.glmer1d, type='diag')
## $Nest
Conclusions:
performance::check_model(owls.glmer1d)
Conclusions:
It is also possible to predict which modelling family would be the most suitable for the data. This can be attempted via an experimental routine in the performance package that uses random forests to classify a range of possible distributions.
performance::check_distribution(owls.glmer1d)
Conclusions:
We can explore over-dispersion more formally:
performance::check_overdispersion(owls.glmer1d)
## # Overdispersion test
##
## dispersion ratio = 4.211
## Pearson's Chi-Squared = 2463.718
## p-value = < 0.001
Conclusions:
Over-dispersion could be caused by numerous factors:
We can explore zero-inflation directly:
performance::check_zeroinflation(owls.glmer1d)
## # Check for zero-inflation
##
## Observed zeros: 156
## Predicted zeros: 44
## Ratio: 0.28
Conclusions:
owls.resid = simulateResiduals(owls.glmer1d, plot=TRUE, integerResponse = TRUE)
Conclusions:
Perhaps we should specifically explore zero-inflation.
testZeroInflation(owls.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 2.6233, p-value < 2.2e-16
## alternative hypothesis: two.sided
Conclusions:
The data were collected at various times throughout the night. It is possible that this could lead to patterns of dependency that are not already accounted for. For example, perhaps observations that are collected at similar time of the night (within a given nest) have more similar residuals than those at very different time of the night. We can explore whether there are any temporal autocorrelation patterns.
testTemporalAutocorrelation(owls.resid, time=owls$ArrivalTime)
## Error in testTemporalAutocorrelation(owls.resid, time = owls$ArrivalTime): testing for temporal autocorrelation requires unique time values - if you have several observations per time value, either use the recalculateResiduals function to aggregate residuals per time step, or extract the residuals from the fitted object, and plot / test each of them independently for temporally repeated subgroups (typical choices would be location / subject etc.). Note that the latter must be done by hand, outside testSpatialAutocorrelation.
owls.resid1 <- recalculateResiduals(owls.resid, group=interaction(owls$ArrivalTime, owls$Nest), aggregateBy = mean)
testTemporalAutocorrelation(owls.resid1, time=unique(owls$ArrivalTime))
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 1.997, p-value = 0.9787
## alternative hypothesis: true autocorrelation is not 0
Conclusions:
plot_model(owls.glmmTMB1d, type='diag')
## $Nest
Conclusions:
performance::check_model(owls.glmmTMB1d)
Conclusions:
It is also possible to predict which modelling family would be the most suitable for the data. This can be attempted via an experimental routine in the performance package that uses random forests to classify a range of possible distributions.
performance::check_distribution(owls.glmmTMB1d)
Conclusions:
We can explore over-dispersion more formally:
performance::check_overdispersion(owls.glmmTMB1d)
## # Overdispersion test
##
## dispersion ratio = 4.178
## Pearson's Chi-Squared = 2460.744
## p-value = < 0.001
Conclusions:
Over-dispersion could be caused by numerous factors:
We can explore zero-inflation directly:
performance::check_zeroinflation(owls.glmmTMB1d)
## # Check for zero-inflation
##
## Observed zeros: 156
## Predicted zeros: 44
## Ratio: 0.28
Conclusions:
owls.resid = simulateResiduals(owls.glmmTMB1d, plot=TRUE, integerResponse = TRUE)
Conclusions:
Perhaps we should specifically explore zero-inflation.
testZeroInflation(owls.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 2.6152, p-value < 2.2e-16
## alternative hypothesis: two.sided
Conclusions:
The data were collected at various times throughout the night. It is possible that this could lead to patterns of dependency that are not already accounted for. For example, perhaps observations that are collected at similar time of the night (within a given nest) have more similar residuals than those at very different time of the night. We can explore whether there are any temporal autocorrelation patterns.
testTemporalAutocorrelation(owls.resid, time=owls$ArrivalTime)
## Error in testTemporalAutocorrelation(owls.resid, time = owls$ArrivalTime): testing for temporal autocorrelation requires unique time values - if you have several observations per time value, either use the recalculateResiduals function to aggregate residuals per time step, or extract the residuals from the fitted object, and plot / test each of them independently for temporally repeated subgroups (typical choices would be location / subject etc.). Note that the latter must be done by hand, outside testSpatialAutocorrelation.
owls.resid1 <- recalculateResiduals(owls.resid, group=interaction(owls$ArrivalTime, owls$Nest), aggregateBy = mean)
testTemporalAutocorrelation(owls.resid1, time=unique(owls$ArrivalTime))
##
## Durbin-Watson test
##
## data: simulationOutput$scaledResiduals ~ 1
## DW = 2.0085, p-value = 0.9397
## alternative hypothesis: true autocorrelation is not 0
Conclusions:
Conclusions:
glmer(), so we will proceed with glmmTMB() only.Zero-inflated vs hurdle models
owls.glmmTMB3 <- glmmTMB(NCalls ~ FoodTreatment*SexParent + offset(log(BroodSize)) +
(FoodTreatment*SexParent|Nest),
ziformula=~1, data=owls, family=poisson(link='log'),
REML=TRUE)
#OR
owls.glmmTMB3 <- update(owls.glmmTMB1d, ziformula=~1)
Model validation
owls.resid = simulateResiduals(owls.glmmTMB3, plot=TRUE, integerResponse = TRUE)
testZeroInflation(owls.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0055, p-value = 0.96
## alternative hypothesis: two.sided
testDispersion(owls.resid)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.44372, p-value = 0.112
## alternative hypothesis: two.sided
testUniformity(owls.resid)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.078674, p-value = 0.001204
## alternative hypothesis: two-sided
testQuantiles(owls.resid)
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.0002957
## alternative hypothesis: both
testResiduals(owls.resid)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.078674, p-value = 0.001204
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.44372, p-value = 0.112
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 3, observations = 599, p-value = 0.96
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.02003339
## sample estimates:
## outlier frequency (expected: 0.006110183639399 )
## 0.005008347
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.078674, p-value = 0.001204
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.44372, p-value = 0.112
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 3, observations = 599, p-value = 0.96
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.02003339
## sample estimates:
## outlier frequency (expected: 0.006110183639399 )
## 0.005008347
owls.glmmTMB4 <- glmmTMB(NCalls ~ FoodTreatment*SexParent + offset(log(BroodSize)) +
(FoodTreatment*SexParent|Nest),
ziformula=~FoodTreatment*SexParent, data=owls,
family=poisson(link='log'), REML=TRUE)
Model validation
owls.resid = simulateResiduals(owls.glmmTMB4, plot=TRUE, integerResponse = TRUE)
testZeroInflation(owls.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0004, p-value = 1
## alternative hypothesis: two.sided
testDispersion(owls.resid)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.65935, p-value = 0.264
## alternative hypothesis: two.sided
testUniformity(owls.resid)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.064467, p-value = 0.01376
## alternative hypothesis: two-sided
testQuantiles(owls.resid)
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.01343
## alternative hypothesis: both
testResiduals(owls.resid)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.064467, p-value = 0.01376
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.65935, p-value = 0.264
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 599, p-value = 0.6
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.03347245
## sample estimates:
## outlier frequency (expected: 0.00686143572621035 )
## 0.001669449
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.064467, p-value = 0.01376
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.65935, p-value = 0.264
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 1, observations = 599, p-value = 0.6
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.03347245
## sample estimates:
## outlier frequency (expected: 0.00686143572621035 )
## 0.001669449
owls.glmmTMB5 <- glmmTMB(NCalls ~ FoodTreatment*SexParent + offset(log(BroodSize)) +
(FoodTreatment*SexParent|Nest),
ziformula = ~1, data=owls,
family=nbinom2(link='log'), REML=TRUE)
## Warning in fitTMB(TMBStruc): Model convergence problem; non-positive-definite
## Hessian matrix. See vignette('troubleshooting')
## Warning in fitTMB(TMBStruc): Model convergence problem; false convergence (8).
## See vignette('troubleshooting')
Model validation
owls.resid = simulateResiduals(owls.glmmTMB5, plot=TRUE, integerResponse = TRUE)
testZeroInflation(owls.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0046, p-value = 0.984
## alternative hypothesis: two.sided
testDispersion(owls.resid)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.47439, p-value = 0.008
## alternative hypothesis: two.sided
testUniformity(owls.resid)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.051742, p-value = 0.08093
## alternative hypothesis: two-sided
testQuantiles(owls.resid)
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.0002411
## alternative hypothesis: both
testResiduals(owls.resid)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.051742, p-value = 0.08093
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.47439, p-value = 0.008
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.02090985
## sample estimates:
## outlier frequency (expected: 0.00649415692821369 )
## 0
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.051742, p-value = 0.08093
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.47439, p-value = 0.008
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.1
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.02090985
## sample estimates:
## outlier frequency (expected: 0.00649415692821369 )
## 0
owls.glmmTMB6 <- glmmTMB(NCalls ~ FoodTreatment*SexParent + offset(log(BroodSize)) +
(FoodTreatment+SexParent|Nest),
ziformula = ~FoodTreatment*SexParent, data=owls,
family=nbinom2(link='log'), REML=TRUE)
Model validation
owls.resid = simulateResiduals(owls.glmmTMB6, plot=TRUE, integerResponse = TRUE)
testZeroInflation(owls.resid)
##
## DHARMa zero-inflation test via comparison to expected zeros with
## simulation under H0 = fitted model
##
## data: simulationOutput
## ratioObsSim = 1.0046, p-value = 1
## alternative hypothesis: two.sided
testDispersion(owls.resid)
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.57518, p-value = 0.048
## alternative hypothesis: two.sided
testUniformity(owls.resid)
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.060199, p-value = 0.02603
## alternative hypothesis: two-sided
testQuantiles(owls.resid)
##
## Test for location of quantiles via qgam
##
## data: simulationOutput
## p-value = 0.03123
## alternative hypothesis: both
testResiduals(owls.resid)
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.060199, p-value = 0.02603
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.57518, p-value = 0.048
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.22
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01510851
## sample estimates:
## outlier frequency (expected: 0.00582637729549249 )
## 0
## $uniformity
##
## One-sample Kolmogorov-Smirnov test
##
## data: simulationOutput$scaledResiduals
## D = 0.060199, p-value = 0.02603
## alternative hypothesis: two-sided
##
##
## $dispersion
##
## DHARMa nonparametric dispersion test via sd of residuals fitted vs.
## simulated
##
## data: simulationOutput
## ratioObsSim = 0.57518, p-value = 0.048
## alternative hypothesis: two.sided
##
##
## $outliers
##
## DHARMa bootstrapped outlier test
##
## data: simulationOutput
## outliers at both margin(s) = 0, observations = 599, p-value = 0.22
## alternative hypothesis: two.sided
## percent confidence interval:
## 0.00000000 0.01510851
## sample estimates:
## outlier frequency (expected: 0.00582637729549249 )
## 0
AICc(owls.glmmTMB3, owls.glmmTMB4, owls.glmmTMB5, owls.glmmTMB6)
Conclusions:
owls.glmmTMB4).glmmTMB6 as an additional exercise.plot_model(owls.glmmTMB4, type='eff', terms=c('FoodTreatment', 'SexParent'))
These predictions appear to be based on the mean BroodSize of approximately 4.39.
plot(allEffects(owls.glmmTMB4), multiline=TRUE, ci.style='bars')
These predictions also appear to be based on the mean BroodSize, although the documentation seems to suggest that allEffects() might not deal with the offsets the way we have used them (as a function in the formula) correctly.
ggpredict(owls.glmmTMB4, terms=c('FoodTreatment', 'SexParent')) %>% plot
This seems to deal with the offset incorrectly. For the purpose of prediction, the offset seems to be set at the value of the first BroodSize (on the response scale). This is incorrect for two reasons:
ggemmeans() can accommodate the offset correctly. There are two sensible choices:
#off<-owls %>% group_by(SexParent, FoodTreatment) %>% summarize(Mean=mean(BroodSize))
off<-owls %>% summarize(Mean=mean(BroodSize))
as.numeric(off)
## [1] 4.392321
ggemmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, offset=log(off$Mean)) %>% plot
ggemmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, offset=0) %>% plot
summary(owls.glmmTMB4)
## Family: poisson ( log )
## Formula:
## NCalls ~ FoodTreatment * SexParent + offset(log(BroodSize)) +
## (FoodTreatment * SexParent | Nest)
## Zero inflation: ~FoodTreatment * SexParent
## Data: owls
##
## AIC BIC logLik deviance df.resid
## 3849.1 3928.3 -1906.6 3813.1 585
##
## Random effects:
##
## Conditional model:
## Groups Name Variance Std.Dev. Corr
## Nest (Intercept) 0.1618 0.4022
## FoodTreatmentSatiated 1.4356 1.1982 0.14
## SexParentMale 0.1035 0.3216 -0.44 -0.28
## FoodTreatmentSatiated:SexParentMale 0.6257 0.7910 0.20 -0.36 -0.67
##
##
##
##
##
## Number of obs: 599, groups: Nest, 27
##
## Conditional model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) 0.86455 0.08960 9.649 < 2e-16 ***
## FoodTreatmentSatiated -0.82476 0.27748 -2.972 0.00296 **
## SexParentMale -0.12590 0.08708 -1.446 0.14824
## FoodTreatmentSatiated:SexParentMale 0.22165 0.23233 0.954 0.34007
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Zero-inflation model:
## Estimate Std. Error z value Pr(>|z|)
## (Intercept) -1.3762 0.2260 -6.090 1.13e-09 ***
## FoodTreatmentSatiated 0.6988 0.3197 2.186 0.0288 *
## SexParentMale -0.8363 0.3331 -2.511 0.0121 *
## FoodTreatmentSatiated:SexParentMale 0.6956 0.4471 1.556 0.1198
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Conclusions:
tidy(owls.glmmTMB4, conf.int=TRUE)
## or on the response scale
tidy(owls.glmmTMB4, conf.int=TRUE, exponentiate = TRUE)
tidy(owls.glmmTMB4, conf.int=TRUE, exponentiate = TRUE) %>% kable
| effect | component | group | term | estimate | std.error | statistic | p.value | conf.low | conf.high |
|---|---|---|---|---|---|---|---|---|---|
| fixed | cond | NA | (Intercept) | 2.3739302 | 0.2127051 | 9.648917 | 0.0000000 | 1.9915900 | 2.8296710 |
| fixed | cond | NA | FoodTreatmentSatiated | 0.4383406 | 0.1216319 | -2.972291 | 0.0029559 | 0.2544592 | 0.7551014 |
| fixed | cond | NA | SexParentMale | 0.8817058 | 0.0767773 | -1.445792 | 0.1482355 | 0.7433660 | 1.0457905 |
| fixed | cond | NA | FoodTreatmentSatiated:SexParentMale | 1.2481300 | 0.2899739 | 0.954029 | 0.3400690 | 0.7915949 | 1.9679618 |
| fixed | zi | NA | (Intercept) | 0.2525244 | 0.0570675 | -6.089913 | 0.0000000 | 1.3128496 | 1.8116620 |
| fixed | zi | NA | FoodTreatmentSatiated | 2.0113883 | 0.6430705 | 2.185777 | 0.0288319 | 2.1274547 | 6.6973290 |
| fixed | zi | NA | SexParentMale | 0.4332912 | 0.1443374 | -2.510653 | 0.0120508 | 1.1844919 | 1.8423425 |
| fixed | zi | NA | FoodTreatmentSatiated:SexParentMale | 2.0048423 | 0.8963908 | 1.555682 | 0.1197838 | 1.5392138 | 4.2663494 |
| ran_pars | cond | Nest | sd__(Intercept) | 0.4021858 | NA | NA | NA | -0.3646182 | 0.5639225 |
| ran_pars | cond | Nest | sd__FoodTreatmentSatiated | 1.1981861 | NA | NA | NA | -0.6358481 | 0.1582352 |
| ran_pars | cond | Nest | sd__SexParentMale | 0.3216476 | NA | NA | NA | -0.2441513 | 0.7805384 |
| ran_pars | cond | Nest | sd__FoodTreatmentSatiated:SexParentMale | 0.7909938 | NA | NA | NA | -0.2888018 | 0.4834026 |
| ran_pars | cond | Nest | cor__(Intercept).FoodTreatmentSatiated | 0.1441205 | NA | NA | NA | -0.3664646 | 0.7170364 |
| ran_pars | cond | Nest | cor__(Intercept).SexParentMale | -0.4407330 | NA | NA | NA | -0.0053304 | 0.6138122 |
| ran_pars | cond | Nest | cor__(Intercept).FoodTreatmentSatiated:SexParentMale | 0.1983013 | NA | NA | NA | -1.8191758 | -0.9333190 |
| ran_pars | cond | Nest | cor__FoodTreatmentSatiated.SexParentMale | -0.2784433 | NA | NA | NA | 0.0721958 | 1.3254545 |
| ran_pars | cond | Nest | cor__FoodTreatmentSatiated.FoodTreatmentSatiated:SexParentMale | -0.3591540 | NA | NA | NA | -1.4892457 | -0.1834447 |
| ran_pars | cond | Nest | cor__SexParentMale.FoodTreatmentSatiated:SexParentMale | -0.6705776 | NA | NA | NA | -0.1807597 | 1.5718905 |
Conclusions:
# warning this is only appropriate for html output
sjPlot::tab_model(owls.glmmTMB4, show.se=TRUE, show.aic=TRUE)
| N Calls | |||||
|---|---|---|---|---|---|
| Predictors | Incidence Rate Ratios | std. Error | CI | p | |
| Count Model | |||||
| (Intercept) | 2.37 | 0.09 | 1.99 – 2.83 | <0.001 | |
| FoodTreatment [Satiated] | 0.44 | 0.28 | 0.25 – 0.76 | 0.003 | |
| SexParent [Male] | 0.88 | 0.09 | 0.74 – 1.05 | 0.148 | |
|
FoodTreatment [Satiated] * SexParent [Male] |
1.25 | 0.23 | 0.79 – 1.97 | 0.340 | |
| Zero-Inflated Model | |||||
| (Intercept) | 0.25 | 0.23 | 0.16 – 0.39 | <0.001 | |
| FoodTreatment [Satiated] | 2.01 | 0.32 | 1.07 – 3.76 | 0.029 | |
| SexParent [Male] | 0.43 | 0.33 | 0.23 – 0.83 | 0.012 | |
|
FoodTreatment [Satiated] * SexParent [Male] |
2.00 | 0.45 | 0.83 – 4.82 | 0.120 | |
| Random Effects | |||||
| σ2 | 0.22 | ||||
| τ00 Nest | 0.16 | ||||
| τ11 Nest.FoodTreatmentSatiated | 1.44 | ||||
| τ11 Nest.SexParentMale | 0.10 | ||||
| τ11 Nest.FoodTreatmentSatiated:SexParentMale | 0.63 | ||||
| ρ01 | 0.14 | ||||
| -0.44 | |||||
| 0.20 | |||||
| ICC | 0.77 | ||||
| N Nest | 27 | ||||
| Observations | 599 | ||||
| Marginal R2 / Conditional R2 | 0.111 / 0.799 | ||||
| AIC | 3849.144 | ||||
options(width=100)
tidy(owls.glmmTMB8, conf.int=TRUE)
plogis(-1.3705)
exp(-1.3705)
tidy(owls.glmmTMB8, effects='fixed', conf.int=TRUE, exponentiate=TRUE)
<class=‘HIDDEN’>
#r.squaredGLMM(owls.glmmTMB4)
performance::r2_nakagawa(owls.glmmTMB4)
## # R2 for Mixed Models
##
## Conditional R2: 0.799
## Marginal R2: 0.111
owls.grid = with(owls, list(FoodTreatment=levels(FoodTreatment),
SexParent=levels(SexParent)))
newdata = emmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, at=owls.grid,
offset=0, type='response') %>%
as.data.frame
head(newdata)
ggplot(newdata, aes(y=rate, x=FoodTreatment)) +
geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL, color=SexParent),
position=position_dodge(width=0.2)) +
scale_y_continuous('Number of sibling negotiations per chick') +
theme_bw()
##OR if we want to express this for the average brood size
owls.grid = with(owls, list(FoodTreatment=levels(FoodTreatment),
SexParent=levels(SexParent)))
newdata = emmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, at=owls.grid,
offset=log(mean(owls$BroodSize)), type='response') %>%
as.data.frame
head(newdata)
ggplot(newdata, aes(y=rate, x=FoodTreatment)) +
geom_pointrange(aes(ymin=lower.CL, ymax=upper.CL, color=SexParent),
position=position_dodge(width=0.2)) +
scale_y_continuous('Number of sibling negotiations per nest') +
theme_bw()
newdata=tidy(owls.glmmTMB3, effects='fixed', conf.int=TRUE, exponentiate=TRUE) %>%
mutate(Model='zip (simple zi)') %>%
bind_rows(
tidy(owls.glmmTMB4, effects='fixed', conf.int=TRUE, exponentiate=TRUE) %>%
mutate(Model='zip (complex zi)')
) %>%
bind_rows(
tidy(owls.glmmTMB5, effects='fixed', conf.int=TRUE, exponentiate=TRUE) %>%
mutate(Model='zinb (simple zi)')
) %>%
bind_rows(
tidy(owls.glmmTMB6, effects='fixed', conf.int=TRUE, exponentiate=TRUE) %>%
mutate(Model='zinb (complex zi)')
) %>%
mutate(Model=factor(Model, levels=c('zip (simple zi)', 'zip (complex zi)',
'zinb (simple zi)', 'zinb (complex zi)')),
Cond=interaction(component, term)) %>%
arrange(component, term) %>%
mutate(Cond=factor(Cond, levels=rev(unique(Cond))))
## Warning in sqrt(diag(vv)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
## Warning in sqrt(diag(vv)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
## Warning in sqrt(diag(object$cov.fixed)): NaNs produced
ggplot(newdata, aes(y=estimate, x=Cond, color=Model)) +
geom_pointrange(aes(ymin=conf.low, ymax=conf.high), position=position_dodge(width=0.2)) +
coord_flip()
newdata = emmeans(owls.glmmTMB3, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
as.data.frame %>% mutate(Model='zip (simple zi)', response=rate) %>%
bind_rows(
emmeans(owls.glmmTMB4, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
as.data.frame %>% mutate(Model='zip (complex zi)', response=rate)
) %>%
bind_rows(
emmeans(owls.glmmTMB5, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
as.data.frame %>% mutate(Model='zinb (simple zi)', response=response)
) %>%
bind_rows(
emmeans(owls.glmmTMB6, ~FoodTreatment+SexParent, offset=0, at=owls.grid, type='response') %>%
as.data.frame %>% mutate(Model='zinb (complex zi)', response=response)
) %>%
mutate(Model=factor(Model, levels=c('zip (simple zi)', 'zip (complex zi)',
'zinb (simple zi)', 'zinb (complex zi)')))
head(newdata)
ggplot(newdata, aes(y=response, x=FoodTreatment)) +
geom_pointrange(aes(color=SexParent, ymin=lower.CL, ymax=upper.CL),
position=position_dodge(width=0.2)) +
facet_wrap(~Model, nrow=1)